function [best_p_aic, best_model_aic, best_p_bic, best_model_bic, ...
    aic_values, bic_values] = choose_best_ar_model(y, p_max)
    % choose_best_ar_model: Chooses the best AR(p) model based on both AIC and BIC.
    %
    % Inputs:
    %   - y: A vector representing the time series {y0, y1, y2, ..., yt}
    %   - p_max: The maximum AR order (p) to consider
    %
    % Outputs:
    %   - best_p_aic: The order of the best AR(p) model selected according to AIC
    %   - best_model_aic: The estimated AR model object (arima model) with the best AIC
    %   - best_p_bic: The order of the best AR(p) model selected according to BIC
    %   - best_model_bic: The estimated AR model object (arima model) with the best BIC
    %   - aic_values: A vector of AIC values for models from p=0 to p_max
    %   - bic_values: A vector of BIC values for models from p=0 to p_max

    disp_result = false;

    % Initialize matrices to store log-likelihood and model orders
    LogL = zeros(p_max+1, 1);  % Log-likelihood matrix
    PQ   = zeros(p_max+1, 1);  % Model complexity matrix (p)

    % Number of observations in the time series
    T = length(y);

    % Loop over all AR orders (p) from 0 to p_max
    for p = 0:p_max
        try
            % Estimate ARIMA(p,0,0) model (AR(p) model)
            if p == 0
                Mdl = arima('Constant', 0);
            else
                Mdl = arima('ARLags', p, 'Constant', 0);
            end
            [~, ~, LogL(p+1)] = estimate(Mdl, y, 'Display', 'off');  % Log-likelihood of the model  
        catch
            % If model estimation fails (e.g., non-convergence), assign -Inf log-likelihood
            LogL(p+1) = -Inf;
        end
        PQ(p+1) = p; % Store the number of estimated parameters
    end

    % Exclude invalid models (where logL or pq are Inf or -Inf)
    valid_indices = ~isinf(LogL);

    % Check if all models are invalid
    if sum(valid_indices) == 0
        error('All AR(p) models failed to converge. No valid models to choose from.');
    end

    % Filter out invalid log-likelihoods and pq values
    logL_valid = LogL(valid_indices);
    pq_valid = PQ(valid_indices);

    % Calculate AIC and BIC for valid models
    [aic, bic] = aicbic(logL_valid, pq_valid+1, T);

    % Initialize the AIC and BIC value vectors
    aic_values = Inf(p_max+1, 1);  % Initialize with Inf
    bic_values = Inf(p_max+1, 1);  % Initialize with Inf

    % Fill in valid AIC and BIC values
    aic_values(valid_indices) = aic;
    bic_values(valid_indices) = bic;

    % Find the minimum AIC and BIC and their corresponding p values
    [min_aic_value, best_index_aic] = min(aic_values);
    [min_bic_value, best_index_bic] = min(bic_values);

    % Output the best AR orders (p) based on AIC and BIC
    best_p_aic = best_index_aic - 1;  % Adjust for MATLAB's 1-based indexing
    best_p_bic = best_index_bic - 1;

    % Estimate the best models using the selected p values
    % Estimate ARIMA(p,0,0) model (AR(p) model)
    if best_p_aic == 0
        best_model_aic = arima('Constant',0);
    else
        best_model_aic = arima('ARLags', best_p_aic, 'Constant', 0);
    end
    if best_p_bic == 0
        best_model_bic = arima('Constant',0);
    else
        best_model_bic = arima('ARLags', best_p_bic, 'Constant', 0);
    end
    
    % Estimate the best models (this is just to store the model, no output needed)
    best_model_aic = estimate(best_model_aic, y, 'Display', 'off');
    best_model_bic = estimate(best_model_bic, y, 'Display', 'off');

    % Display the results (optional)
    if disp_result == true
        disp('AIC values:');
        disp(aic_values);
        disp(['Minimum AIC: ', num2str(min_aic_value)]);
        disp(['Best AR order (p) according to AIC: ', num2str(best_p_aic)]);

        disp('BIC values:');
        disp(bic_values);
        disp(['Minimum BIC: ', num2str(min_bic_value)]);
        disp(['Best AR order (p) according to BIC: ', num2str(best_p_bic)]);
    end
end